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Abstract 

We consider the modeling of the effect of unresolved scales, for two-dimensional and 
geophysical flows. We first show that the effect of small scales on a coarse-grained field, 
can be approximated at leading order, by the effect of the strain tensor on the gradient 
of the vorticity, which exactly conserves the energy. We show that this approximation 
would lead to unstable numerical code. In order to propose a stable parameterization, 
while taking into account of these dynamical properties, we apply a maximum entropy 
production principle. The parameterization acts as a selective diffusion proportional 
to the mean strain, in the contraction direction, while conserving the energy. We show 
on numerical computation that the obtained anisotropic relaxation equations give an 
important predictability improvement, with respect to Navier-Stokes, Smagorinsky or 
hyperviscous parameterizations. 
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1 Introduction 



Turbulent flows are characterized by the non-linear interaction of a huge number of vari- 
ables. For very high Reynolds' numbers, the current or foreseen computational capabilities 
only permit to describe part of corresponding degrees of freedom. For instance, for geophys- 
ical flows, the molecular viscosity acts at typical scales of order of the millimeter, whereas a 
current typical resolution for meteorological computations is about 30 to 50 km. The issue 
of the parameterization of the effect of unresolved scale on the actually described ones is 
thus of crucial importance. 

The predominance of the Coriolis force and the density stratification give geophysical 
flows a quasi two-dimensional nature. Their organization is actually three dimensional, but 
may be well approximate by layered models [1]. Moreover, at large scale, an inverse energy- 
cascade process (from small to large scale), typical of two-dimensional flows, is observed for 
geophysical flows. In the following, we address the issue of small scale turbulence param- 
eterization in both the frameworks of the two-dimensional incompressible Euler equation, 
or of the Quasi-Geostrophic model on the other hand. 

Any numerical resolution of the flow dynamics implicitly assumes that an averaged 
value is described. The average of the nonlinear terms should then be expressed in terms 
of the known averages of the single quantities. This is the classical closure problem for flow 
equations [2]. In order to clearly specify the problem, one has to precise the meaning of 
the average. When a statistically stationary situation is identified, a very natural choice 
is an ensemble average on the stationary distribution. For instance, this would be the 
case, for an inverse cascade regime (a small scale energy injection and a large scale drag 
force) where a natural ensemble is given by a stationary probability distribution for this 
process. The proximity of such a stationary situation is implicitly assumed, together with 
other hypothesis, in classical parameterizations such as the Eddy Damped Quasi Markovian 
Model (EDQNM) [3, 4, 5] or the K - e model [4]. An other approach is to specify the 
stochastic properties of the unresolved scales, as for instance in the Direct Interaction 
Approximation (DIA) [7]. 

However, geophysical flows are characterized by the existence of large scales structures 
(jets, cyclone and anticyclones), which break the self- similarity at the base of the description 
of the inverse cascade process. The typical time scales for forcing and dissipation are much 
larger than the inertial ones. In order to study such flows, we will consider two dimensional 
inviscid models in a freely decaying-turbulence regime. In such a case, for large times, 
the flow self-organize in coherent structures. The final state may be described by the 
statistical mechanics of the vorticity (or the potential vorticity) [9, 10]. For such a freely 
decaying turbulence, and in a situation close to the equilibrium or to a quasi-stationary 
state of the inviscid equations, it would be very natural to consider an ensemble average, 
compatible with the observed large scale structure and vorticity distribution. Such an 
ensemble average is for instance implicitly assume in the kinetic approach of Chavanis [11]. 
On the contrary, in this work, we will be interested in situations far from equilibrium or 
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from stationary solutions (vortex merging, filamentation). In such a case, there is no natural 
statistical ensemble to consider. Moreover, we will show that for a given flow evolution, the 
effect of unresolved scales will be dominated by structures (filaments) which are strongly 
correlated with the resolved scales. Given this situation, we will consider the evolution of 
locally averaged quantities, on the grid scale I (average by homogenization by opposition 
of ensemble average) . 

In section 2 we recall the equation of motion and give the formal properties of an 
average by homogenization. We then show that the averaged equations will be determined 
by the knowledge of a turbulent current. In section 3, we give an analytic expression of 
this turbulent current, at leading order in an asymptotic expansion in power of the ratio of 
grid scale on the domain scale l/L. The turbulent current is then expressed in terms of the 
averaged velocity strain tensor. This result has been described in [14] and independently 
in [16, 15]. We discuss the validity of this deterministic approximation. Using numerical 
integration, we show that this expression is actually a very good approximation of the 
actual turbulent current, in situations of decaying turbulence, far from equilibrium. This 
deterministic approximation exactly conserves the energy of the mean flow. 

As the strain tensor is of zero trace, its two eigenvalues are opposite. They respectively 
correspond to a positive and negative viscosity. For this reason, this deterministic approx- 
imation is not suited to model the effect of small scales in a numerical code: it would lead 
to instability of the scheme. In section 4, we obtain a stable algorithm by determining the 
turbulent current which maximize the production of the mixing entropy, while respecting 
both the value of the deterministic current in the stable direction of the strain tensor, and 
the conservation of energy. Such a maximum entropy production principle has already been 
used by Robert and Sommeria [12], but without taking into account for the anisotropy of 
the turbulent current. With this new specification, the obtained equations respect all the 
usable properties of the deterministic current, and have the further property to converge 
towards equilibrium of the statistical mechanics. In section 5, we compare the performance 
of these anisotropic relaxation equations with other models, commonly used to parame- 
terize small scale turbulence, for two-dimensional or geophysical flows: the eddy viscosity 
model, the Smagorinsky model [13], and an hyperviscosity parameterization. In section 5, 
the anisotropic relaxation equation are shown to give the best results, when compared with 
high resolution simulation, both for the Euler and for the Quasi-Geostrophic equations. 

In section 6, we discuss these results and their interest for geophysical flow applications. 

2 The averaged equations of motion 

In this section, we recall the two-dimensional incompressible Euler equations and the Quasi- 
Geostrophic model. We specify the average operator for the resolved dynamics and we give 
the main properties of the averaged equations of motion. 
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The incompressible Euler equations describe the evolution of the velocity u of an incom- 
pressible inviscid fluid. For 2D flows, they are equivalent to the equation for the evolution 
of the vorticity (w = (VA u) .e z ): 



— + u.Vlj = 
at 

uj = —Aip 



(1) 



u = -e z A VY> 



(2) 
(3) 



The vorticity to is advected by the divergence-less velocity u, ip is the stream-function. 
For sake of simplicity, we will consider these equation on a torus D (all the quantities are 
periodic in both variables, with zero global average). All our discussion is easily generalizable 
for any bounded domain D (an impermeability condition u.n, where n is normal to the 
boundary of D would then specify the boundary conditions). 
The Euler equations conserve the energy : 



and all the functionals of the form Cf(u>) = j D dxf{u) for any function /. The Quasi- 
Geostrophic (QG) equation, which is the simplest model of geophysical flows [1], are : 



We note the formal similarity with the Euler equation. The potential vorticity q is then 
advected. The determination of the stream function, is then given by the inversion of the 
equation (6), and the velocity is still given by (3). The energy for the QG dynamics is 
E = 1/2 J D dr [u 2 + t(j 2 /R 2 ] = 1/2 J D drqip. One clearly see the strong analogies between 
the QG and the Euler equations. From now on, we will deal only with the Euler equation. 
We will precise when the generalization for the QG equation is not straightforward. 

Our aim is to obtain the equation describing the evolution of the spatial average of the 
velocity and of the vorticity. We define the average / of any quantity / by : 



This operator is the convolution with G. We suppose that / drG(r) = 1 and that G 
is invariant under any rotation (G(r) only depends on the modulus of r.) The qualitative 
properties assumed for G are to regularize the field at a fixed scale I. For instance, in 
section 3, we will use in the 27r-periodic domain the kernel G l whose Fourier components 




(4) 




(5) 



(6) 
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are = ^ exp ^— , k G Z 2 (when Z is much smaller than 27r, G l (r) is close to 

2^pj exp (—^p^j)- The average operator verifies the linearity property : / + Xg = f + Xg 
and the commutation with the differential operators : dif = dif. We note that it does not 
verify the properties f = f and fg = fg as would be the case for an ensemble average. In 
section 3, we will show that the dominant contribution of the turbulent current comes from 
the non- verification of these two properties. 

By averaging (1), we obtain the averaged Euler equation : 

— + u.Vcj = — V. J w with J w = lou — iou (7) 

uJ = -Alp (8) 

u = —e z A VV' (9) 

We note that the correlation term is the divergence of a current J w , that we will call 
the turbulent current. We note that the average Euler equations have the same symmetry 
property than the initial Euler equations. 

The aim of our study is to determine an expression of the turbulent current J w in terms 
of the average quantities. This closure problem is equivalent to the usual one for the velocity 
correlations, when the evolution of the velocity is considered. 



3 Anisotropy of the turbulent current 

In this section, we first derive an approximation of the turbulent current, as the leading 
order of an expansion in power of the homogenization scale I (the actual expansion param- 
eter is l/L, where L is a typical size of the domain D). We discuss the main properties 
of this deterministic expression. The turbulent current is anisotropic. It is related to the 
symmetric part of the strain tensor for the averaged velocity. It exactly conserves the en- 
ergy of the averaged fields. Using typical vorticity fields of freely decaying turbulence, we 
numerically compute the actual turbulent current, in order to study the validity of this first 
order approximation. 



3.1 A leading order deterministic approximation 

In order to obtain an approximation of the turbulent current J w , we decompose the 
vorticity field as the sum of its average U and of fluctuations Q, to being defined by uj = uj+uj. 
Similarly u = u + u. We then obtain for the turbulent current J w (7) : 

J u e^-uu= (uU-uu) + wH + wu + wu (10) 

" l ' 2 3 4 

1 * .. ' 
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We note that the first term (double average), the second term (correlation between the 
vorticity fluctuations and the mean velocity) and the third term (correlation between the 
velocity fluctuations and the mean velocity) are all identically equal to zero for an ensemble 
average. The decomposition (10) is the commonly used one. However, for reason that shall 
become clear we will use : 

= Jd + Js with 3d = am + dm — uJu et 3 S = am + am (11) 

The first term 3d, which does not depend on the velocity fluctuations, will be called the 
deterministic component of the turbulent current (we will show that it can be expressed, to 
a good approximation, in terms of the averaged quantities). The second term depends on 
the velocity fluctuations. It will be called the stochastic component of the turbulent current. 
As the order of magnitude of u is I times the order of magnitude of u, we anticipate that 
the stochastic current will be negligible at leading order in I. 

In order to relate the deterministic part of the turbulent current 3d = am — am (11) to 
the strain tensor for the averaged velocity, we expand u in Taylor series around a point O. By 
noting that the zeroth order terms are null, we obtain at the leading order in I : 3d (O) = 
Ux (O) d x ll^ (O) + uJy (O) dyU^ (O) . Using that the the average operator in equivalent to 
a Gaussian for small I, we note that a part-integration leads to ZJx (O) = l 2 d x uj(0) and 
uJy (O) = l 2 d y uj{0). Using this result, we obtain that the deterministic part of the turbulent 
current is equal, at leading order, to the action of the strain tensor on the vorticity gradient : 

3d = / 2 SVoJ (order 1) with E = ( ^ ) 

As the typical variation length for the mean vorticity uJ is I, its gradient is of order 1/7, and 
this expression is actually of order 1. This expansion may be led to any order in I. At each 
order it involves derivatives of uJ at larger orders. In section 3.2 we will study numerically 
this approximation for vorticity fields typical of freely decaying turbulence, in the regime 
of vortex merging and filamentation. 

The antisymmetric part of the strain tensor, X! an (, characterizes the local rotation of 

the fluid. It is linked to the vorticity by : T: ant = 2 ( ^_ ^ j . It is easily verified that 

V. (S ant Vo;) = 0. Thus, whereas it is present in the leading order approximation for J^, 
this antisymmetric part has no effect on the evolution of the vorticity. We have : 

V.J^V.PV.V*) (order!) with ^ = ( (a /2 ^ ^ ) 

(12) 

The symmetric part of the strain tensor may be diagonalized with orthogonal eigenvectors. 
As the fluid is incompressible, it has a null trace and thus has two opposite eigenvalues a 
et —a, with eigenvectors e CT and e_o- respectively, a is the mean strain : 
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a = ^/-det (Ss-ym) = \J-d x u x d y Uy + (d y u x - d x u y ) 2 /4 



Let us consider the operator u — ► J D dr uJV. (5] sym Vu;). The contraction direction of 
the strain tensor e CT corresponds to the definite negative part of this operator : for any 
uJ, J D dr ZuV. (e a e a .VZJ) < (by part integration). It thus acts as a positive viscosity. 
Conversely, the stretching direction of the strain e_ CT tensor will correspond to the positive 
part of this operator. It thus act as a negative viscosity We will note ^ S ym = + ^f ym 
with ^f ym = <Je a T e CT and ^f ym = —cre_ a T e_ CT where T is the transposition. From a 
numerical point of view, the negative part of the symmetric part of the strain tensor ^f ym 
will have a stabilization effect whereas the positive part 5]^ m will have a destabilization 
effect. In a numerical code, this last term will lead to instabilities. For this reason, the 
approximate expression (12) does not suit for a direct parameterization of the effect of 
small scales on large scales. 

To conclude this analysis, we note that the leading order approximation of the tur- 
bulent current (12) exactly conserves the energy of the averaged field : E = \ J D dru 2 = 
\ J D drujip. The conservation of the energy is a key property of two-dimensional turbulence. 
In section 3.2, we will show numerically that the actual turbulent current (not the leading 
order approximation) is responsible for a very small flux of energy towards the large scale 
(E slightly increase). This is in accordance with the phenomenology of two-dimensional 
freely decaying turbulence. 

We have described a leading order approximation for the turbulent current. This approx- 
imation is deterministic (it can be expressed explicitly in terms of the averaged quantities). 
The turbulent current is the application of the strain tensor on the gradient of the mean 
vorticity. This approximation exactly conserves the energy. In the following section we will 
show numerically that this approximation is a go one for typical vorticity fields. However, 
this analytical expression can not be used to design a stable numerical scheme. In section 
4, we will propose a stable scheme that respects the main properties of this approximation. 

3.2 Numerical analysis of the turbulent current 

For a given vorticity field, it is possible to compute numerically the averaged field at a 
given scale. From it, one can then compute the turbulent current. The aim of this section 
is to describe qualitatively the turbulent current, and to study the validity of the leading 
order approximation, for vorticity fields which are typical of freely decaying turbulence. 

In order to obtain freely decaying turbulence typical vorticity fields, we use results of 
numerical simulations. For all numerical computations of this article, we will use a pseudo- 
spectral code, with an order 3 Adam-Bradsforth scheme for the temporal discretisation, on 
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a domain periodic in both direction. The result are presented for a vorticity field obtained 
by numerically solving the Navier Stokes equations : 

^ + u.Vuj = i/Aw (13) 
ot 

For a given resolution, the viscosity v is chosen as small as possible, such that the computa- 
tion is stable. The Navier Stokes equation verify a maximum principle : the maxima (resp. 
minima) of the vorticity must always decrease (resp. increase) in time. In order to test the 
stability of the computation, we check for any computation the extrema of the vorticity. 
We have studied the turbulent current for several situations. We discuss the results for the 
vorticity field represented on Fig. 1. It has been obtained by considering an initial condi- 
tion made of random patches of potential vorticity. The typical vorticity and velocity are of 
order 1. The parameter of the computation are v = 3.14 10" 5 , Re = 2Ttu max /v = 200000, 
dt = 6.14 10 -4 (time step), resolution 1024X1024. The vorticity field is the one correspond- 
ing to t = 15. At this time, the decrease of the total energy is of order 2%, such that 
the Navier-Stokes approximation may be considered a good approximation of the Euler 
equation. A more precise description of the computation is given in Bouchet (2001). 

This vorticity field is characteristic of the regime of strong filamentation, due to the 
nonlinear interaction between the vortex. Fig. 1 shows that coherent vortices are present 
together with filament structures, at all scale of the field. We will compute the averaged vor- 
ticity field, using the averaging operator described in section 2, with I = 2-7r/128. This cor- 
responds to an homogenization at a scale corresponding to a numerical resolution 128X128. 
The corresponding averaged vorticity field is represented on Fig. 1. 

From the averaged vorticity field, we compute the turbulent current — ojvl — uj u. 
Fig. 2 shows the divergence V.Joj of the turbulent current, associated to the vorticity field 
uj. The comparison of this picture with the vorticity one (Fig. 1) clearly illustrates that 
the divergence of the turbulent current is associated to structures with typical scales of the 
order of the homogenization scale I. Most of these structures are filaments that are stretched 
by the flow strain, passing from resolved to unresolved scales. Fig. 3 shows the mean strain 
a, for the averaged velocity field. The strain is a vectorial property. This figure however 
clearly show that areas with strong strain are associated with strong vorticity gradients 
(alignment of the structures with the strain), and that these areas are the ones that give 
strong contributions to the divergence of the turbulent current. This qualitative description 
is in accordance with the result we have obtained for the leading order approximation of 
the turbulent current by the effect of the strain tensor on the gradient of the vorticity field. 

Let us analyze quantitatively the contribution of the various terms composing the tur- 
bulent current. Table 1 shows the divergence of the terms of the classical decomposition 
(10). This illustrates that the two first terms have a contribution much larger than the total 
turbulent current. We thus note that these two first current are strongly anti-correlated. 
Indeed their sum give a contribution much less important that each of its components. 
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Table 1 also shows that the principal contribution of J s is the correlation between vorticity 
fluctuations and velocity fluctuations. In a numerical simulation with resolved scales up 
to the scale I, this will correspond to the contribution of unresolved scales, for which any 
information is loss. A statistical modeling of this component will then be the more natural. 
That's why we call J s , the stochastic component of the turbulent current. The deterministic 
part dominates the stochastic part, by a factor 10 to 15, for the three studied norms. This 
is in accordance with the leading order approximation we have obtained in section 3.1. 

Let us discuss the results obtained by Laval, Dubrulle and Nazarenko [17], for currents 
computed for an average corresponding to a spectral truncation of the vorticity field. They 
have studied the relative contribution of the currents 2, 3 et 4 (see 10). They have shown that 
the advection of the fluctuations by the average velocity (term 2) dominate the terms 3 and 
4. This is in accordance with the result we have obtained with a Gaussian homogenization. 
Using these results, they have justified an approximation neglecting these last two terms. 
They have then proposed a parameterization of small scale effects, in which the second term 
is determined thanks to a particle description of the fluctuations [18]. We will follow another 
approach : we use a deterministic approximation of the sum of the two terms 1+2=3^. 

Table 2 shows the norms of the difference between and its first (12) and second order 
approximation, in the small l/L expansion. We conclude that, for this vorticity field, in 
the regime of vortex merging and strong filamentation, the deterministic part of the tur- 
bulent current can be approximated by the action of the strain tensor on the gradient of 
the vorticity, up to an error of order of 10%. Table 2 also shows the relative contribution 
of the positive and of the negative parts of the strain tensor. The negative part (Ef ym , 
positive eigenvalue) has a larger contribution than the positive ones (5]^ m , negative eigen- 
value). This illustrates the irreversibility of the flow evolution : due to this property, the 
moments of the vorticity will evolved towards values corresponding to the further mixing of 
the vorticity. We however recall that energy loss associated to the negative part is exactly 
compensated by energy gain associated to the positive part. 



4 Anisotropic relaxation equations 

In section 3 we have obtained a deterministic approximation for the turbulent current 
J w . This approximation is the action of the strain tensor on the gradient of the vorticity. It 
exactly conserves the energy. However the positive part of the strain tensor acts similarly to 
a negative viscosity. This would render any numerical simulation unstable. In this section, in 
order to obtain a parameterization as close as possible to the actual turbulent current, but 
with the essential practical constraint to lead to a stable numerical scheme, we will search 
for a turbulent current which have the property of the deterministic one in the contraction 
direction of the strain tensor (stabilization effect), and which exactly conserves the energy. 

There are several ways to impose such constraints to a parameterization. We choose 
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to determine a parameterization which is in accordance with the statistical properties of 
the mixing of the vorticity. Equilibrium statistical mechanics of the two-dimensional Euler 
equations describes the most probable mixing of the vorticity for a given distribution of 
potential vorticity and Energy. We will then choose the turbulent current which maximize 
the entropy production, while respecting the energy conservation and the effect of the stable 
contribution of the strain tensor. 

A similar maximum entropy production principle (MEPP) has been used by Robert 
and Sommeria [12] in order to parameterize the small scale turbulence. In this work, author 
however impose an isotropic constraint on the turbulent current. We make a derivation very 
close to their one, but we impose the value of the turbulent current in the stable direction 
of the strain tensor. We thus explicitly take into account the anisotropic nature of the 
turbulent current. We do not think that such a maximum entropy production principle is 
based on any physical or theoretical ground. The anisotropic relaxation turbulent current 
will indeed be different from the observed one. We however use this principle to obtain 
stable equation as close as possible from the correct physical ones. 

In section 5, we will show that the obtained anisotropic relaxation equations lead 
to better results than classical parameterizations in order to determine the evolution of 
coarse_grained flow. 

Let us first briefly introduce the results of the equilibrium statistical mechanics of the 
vorticity. In order to simplify this discussion, we will consider an initial condition made of 
many vorticity patches, but with only two values u> = a± and lo = 0,2, occupying respectively 
areas A\ and A<i. Generalization of this discussion and of the following results to any initial 
distribution of vorticity is however straightforward. We note that the initial distribution 
(values a\ and 02, and areas A\ and A2) are conserved by the inviscid dynamics. We now 
consider a coarse-graining of the vorticity at a scale I (local-average at the scale I, which will 
be interpreted in our case, as the numerical resolution scale). At this coarse-grained scale 
(macroscopic description), we will describe the vorticity fields by the local probabilities p(r) 
and (1 — p(r)) to have respectively the vorticity values a\ and 02, at the position r. The 
expression of the coarse-grained vorticity is then : ZJ = a\p + 02(1 — p). The main result 
of the equilibrium statistical mechanics, is the evaluation of the probability to observe a 
local probability field p(r), given the vorticity distribution a given energy E. When the 
scale I goes to zero, the logarithm of the probability of two probability fields is given by the 
entropy S : 

S= dv s(p) with s(p) = — (plogp + (1 — p) log(l — p)) (14) 
Jd 

(see [10, 9], and [20, 21] for a justification). The most probable state is then given by the 
maximization of the entropy for a given energy and vorticity distribution. This equilibrium 
is a stationary state of the dynamical equations. The main hypothesis of the equilibrium 
statistical mechanics is that the complex dynamical evolution of the flow will lead to a state 



10 



close to the equilibrium one. 



We would like to describe the coarse grained evolution of the vorticity. On the contrary to 
the equilibrium case, no clear principle permits to describe this out of equilibrium relaxation. 
The Euler equations advect the two levels of vorticity with the velocity u. This conservation 
equation will lead to a transport equation for p, by the mean velocity u on one hand, and 
by a turbulent current J as a consequence of the unresolved scales : 

d t p + u.Vp = -V.J (15) 

Using the link between the averaged vorticity and the probability : ZJ = a\p + 02(1 — p), 
the coarse-grained vorticity verify equation (7) with = (a\ — ci2)J '. Our aim is thus to 
determine J w . 

As explained in the introduction, we look for a turbulent current which has a determined 
behavior in a direction given by the strain tensor and which conserve the energy. As this 
does not specify the turbulent current, in order to obtain an equation compatible with the 
statistical tendency of the equation to increase the entropy, we will look at equation which 
maximize the entropy production rate, given these dynamical constraints. 

Using (15), we first compute the entropy and energy production rates : 

~aT = ~ dr 7 =w= V 16 

at Jd [ai - lo) (lo - a 2 ) 

and — = / drW-Jc (17) 
at Jd 

We note that (ai — lo) (uJ — 02) > 0. This expression shows that the entropy production is 
maximal when the current has a direction opposite to the one of the vorticity gradient. 
We want to impose an anisotropic constraint on the turbulent current J w ; we put : 

Ci(r) (J^^ +C2(r) (J^^ = i (ig) 

where the two vectors ei and e 2 (depending on r), and the values of C\ and C 2 will be 
specified latter. We search for the current which maximize the entropy production (16), 
while verifying the energy conservation (17), the anisotropic constraint (18) and the global 
constraint on J w : J*£j dv J w — 0. We thus consider the critical points of the functional : 

S-^-E / dr a (r)C l (r)(3 w .e t ) 2 + (e z AV). f drJ^, 
at at ., n JD Jd 

where (3, a(r) and — e 2 A V are Lagrange parameters associated to the constraints. The 
first variations of this functional lead to : 

3 W = - v i + («i - w) (w - a 2 ) -e,AV)) ei (19) 

i=l,2 
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with Ui = — 1/ (aCi (ai — uJ) (uJ — 02)). 

Using the energy constraint (17) and the condition on the zero global average for J^, 
we compute the entropy production : 

dS = f dr E (J "' e ' )2 



dt Jd r^ 2 v% (ai - ui) (u - a 2 ) 



Noting that (cti —uj)(uj — a 2 ) > 0, we conclude that if the diffusivities are strictly posi- 
tive, the entropy variation is actually an entropy production. With strictly positive diffusiv- 
ities, we can also conclude, that if the flow converges towards a stationary state, the entropy 
production must vanish and thus the current must be zero. In such a case, from (19) we 
obtain e 2 A V = V (j3 (cti — 02) ip + log (uJ — 02) — log (ai — a;)^ . The constant e 2 A V being 
a gradient must be equal to zero. This equation can then be integrated to. This yields : 



UJ 



-A?=^-ta„h(^ (« + «)) (20) 



This is the equations of the statistical equilibrium, which is a stationary state of the Euler 
equations. We thus conclude, that if the equations for the coarse-grained field, with tur- 
bulent current (19), converge; they converge towards a statistical equilibrium. We will call 
these equation anisotropic relaxation equations. 



Let us determine f3 and V. For this let us denote A = J2i=i 2 Id d r v i (°i — uj) (uj — a 2 ) ei T e 
<T is the transposition operator). A is thus an order two symmetric tensor), b = J2i=i 2 Id d r u i 

A -b \ 
-*b c ) 

is a 3X3 symmetric matrix. We then express the energy conservation (17) and the property 
of the zero average for : J D drJ^ = 0. This yields three equations to determine the three 
variables (3 and V : 



1 X kJ Ulllj UJ. (llll > J..' WUJ.UJ.WXJ. J../ V J. UilVl I • J. J-. J. kJ UXJ.l_l.kJ (JtXX WX W_WX UVVW U J XJ.XXJ.XW 11 X W U WXXkJ WX I , i-r t 

(2d vector) and c = X)i=i,2 Id d r v i ( a i —TD){u — a 2 ) (Vifj.e^ (scalar). M = ^ 




. , f . . . , / T,i=i,2 Id dr (Voj.e;)e^ \ 

We prove in note 1 that, for strictly positive v\ et v 2 , M is invertible, except for a flow with 

x Let us prove that the symmetric matrix M is invertible. For this, let us consider M as defining a 
quadratic form, and let us prove that it is definite positive. Using the expressions for A, b and c, the 
linearity of the integral and the relation : *ye*eiy = (e^.y) 2 , a straightforward computation leads to : 

(*y x) M ^ Y x ^ = J D dr E i= i, 2 ^ ( ai - w) (w - a 2 ) (e;. (y - xVi/>)) 2 . Given that (01 - w) (w - 02) > 0, 

this computation shows that the quadratic form is positive (> 0). Moreover, noting that ei and e 2 are 

orthonormal, we conclude that if v\ and vi are strictly positives, (*y a;) M f ^ J can be zero only if 



(y — xVV>) is null on the whole domain. This would imply that Vi/> be constant on the domain. With 
periodic boundary conditions, this is possible only if tp = 0, that is for a flow without motion. 
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ip = 0, that is except for an identically zero velocity. 

We have obtained relaxation equations, which have the property to increase the entropy 
while conserving energy. We will now specify the vectors ej and the values of the diffusivity 
V{. For this we impose a diffusivity l 2 a in the stable direction of the strain tensor e a and 
zero in the unstable direction of the strain tensor e a . This choice is the one which allow 
to best approximate the deterministic approximation of the turbulent current, obtained in 
section 3.1, while insuring that the diffusivities are positive. We then obtain the following 
parameterization : 

+ u.VU = V. (l 2 cj [(Vu7 + (01 - 57) (57 - a 2 ) W>) .e a ] e CT ) (22) 

with _ _ 

(5 = I D draV^e a V^.e a ^ _ ^ 

J D dra(ai - 57) (57 - a 2 ) (Vip.e^ 

To obtain this equation for (5 (23), we have used (21) and neglect the effect of the Lagrange 
parameter V (see 19). We have thus assumed V = 0. (in section 5, we will use the relaxation 
equation to model flow evolutions, we have done computation with V^O, using (21), but 
without noting any important difference with the case V = 0) 

The main hypothesis of the model (22) is that the local, destabilizing action of the 
positive part of the strain tensor be replaced by the non local non destabilizing drift term : 
(3 (ai — 57) (57 — 02) Vip. This terms insures the energy conservation. 

From a numerical point of view, the implementation of this equation has the same level 
of complexity as the implementation of the Navier-Stokes equation. All what is further 
needed is the diagonalization of a 2X2 symmetric matrix to determine the strain tensor 
directions and the mean strain a. The computation cost of these operations (of order N, 
where N is the number of degrees of freedom) is negligible with respect to the computation 
of the fast Fourier transform needed to evaluate the nonlinear terms in a pseudo-spectral 
code (of order iVlogiV). 

In next section, we will compare this model with usual parameterizations of two- 
dimensional and geophysical flow computations. We will show that they give more precise 
results for the same resolution (or equivalently for the same computational time). This will 
characterize an important predictability gain for the determination of the flow evolution. 

5 Comparison of the anisotropic relaxation equations with 
commonly used models 

The aim of this chapter is to determine the efficiency of the anisotropic relaxation 
equations when compared with other models of small scale turbulence. We consider the 



13 



regime of vortex merging and filamentation. We study two types of flows : firstly the merging 
of four vortices, which leads to a stable dynamics (the position of the final vortex is always 
on the center of the domain) and secondly the merging of 50 vorticity patches which leads 
to successive vortex merging and for which the positions of the vortices is very sensitive to 
the initial parameterization. In order to compare different models, we will use as a reference, 
high resolution direct numerical computations of the flow. 

We will compare the anisotropic relaxation equation with the following parameteriza- 
tions : 

- direct numerical simulations with the Navier Stokes equations (13) (denoted DNS) 

- with hyperviscous parameterization, where the Laplacian operator is replaced by a 
bi-Laplacian : 

^ + u.Vlu = vAAu (24) 
at 

(denoted HV). Such an hyperviscous parameterization is commonly used in two- 
dimensional and geophysical flow computation. It is not based on any physical ground. 
However, the bi-Laplacian has the effect to stabilize the computation whereas it ac- 
tually acts only on very small scales (higher effective Reynolds number). It usually 
dissipates few energy. It leads to some unphysical spurious description of the small 
scale structures, due to the non-respect of the maximum principle of the Navier Stokes 
equation (the maxima of the vorticity may increase in the hyperviscous case) 

- the Smagorinsky model [13] (denoted SMA). The Smagorinsky model uses an isotropic 
viscosity V. (i/Vw) but with a non-constant viscosity. The viscosity is taken propor- 
tional to the mean shear a : v (r) = l 2 a (r). This parameterization takes into account 
of the importance of the strain for the transport to small scales of the flow struc- 
tures. However the diffusion acts indistinctly in all directions (isotropic). This model 
is commonly used in geophysical fluid computation. 

Concerning the anisotropic relaxation equation, we will alternatively use either the equa- 
tions (22) and (23) (denoted RELANI) or the equations (22) with a local computation of 
the Lagrange parameter (3 for the energy conservation (denoted RELANILOC). In this last 
case we will determine (3 with (22), but with an evaluation of the integrals in subdomains 
Di large compared to the resolution scale I, but smaller than the domain scale. For in- 
stance, for low resolution computation with 256X256 degrees of freedom, we will evaluate 
(3 in 12X12 subdomains. We will also consider the isotropic relaxation equations (Robert 
and Sommeria [12]) (denoted REL) 

For all of these models, the diffusivities v, or the parameters I, are always chosen as 
small as possible for the numerical computation to be stable. 

We quantify the error of the low resolution computation, with respect to reference 
high-resolution computation (field denoted by the subscript re /), by evaluating the ve- 
locity relative error : E r = J D dr (u — u re j) 2 / J D dru^ e j or the vorticity relative error : 
Er-u = J D dr (uj — uJref) 2 / Id ^ ruj ref- (these quantities are evaluated in the Fourier space, 
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using a number of components corresponding to the low resolution fields). 

We first consider an initial condition made of four vortices (Fig. 4). Initially four vorticity 
patches with vorticity a\ = 2.942 are in an ambient vorticity a2 = —0.263. The edges 
of the vortex patches are smoothed to avoid numerical problems for initial times. The 
position of each vortex is slightly moved, in order to break the apparent symmetry of the 
initial condition. We will use as reference computation a Navier-Stokes simulation with 
resolution 1024x1024 (DNS. Euler4. 1024). Numerical parameters are given in table 4. The 
corresponding Reynolds number, based on the domain scale, is then 400 000. Fig. 4 shows the 
vorticity field evolution. This typical experiment of vortex merging show the filamentation 
followed by a very slow organization towards a stationary state. The energy loss between 
times T = and T = 100 is 0.7%. We thus consider this computation as close to the inertial 
limit. 

Fig. 5 shows the evolution of the velocity relative error for 4 computations with reso- 
lutions 256X256 and with respective parameterizations : Navier Stokes (DNS.Euler4.256), 
isotropic relaxation equation (constant diffusivities) (REL.Euler4.256), anisotropic relax- 
ation equations (RELANI.Euler4.256) and anisotropic relaxation with local conservation of 
the energy (RELANILOC.Euler4.256). Table 4 gives the numerical parameters for each of 
these computations. Fig. 5 shows that the maximal relative error is 5%. In this case, the 
numerical evolution is thus easily predictable. After a quasi-exponential growth, the error 
saturates and oscillates together with the position the oscillation of the slightly asymmetric 
vortex. The oscillations reflect the oscillation with respect to the one of the vortex in the 
reference computation. The anisotropic relaxation equation with local conservation of the 
energy give clearly the best results. 

Fig. 6 compares the evolution of the vorticity relative error for the anisotropic re- 
laxation equations with local energy conservation (RELANILOC.Euler4.256), hyperviscos- 
ity (HV.Euler4.256), and the Smagorinsky model (SMA.Euler4.256). For short times, the 
anisotropic relaxations give a slightly better result than the hyperviscous computation. 
However for larger times, the two model are not better one than the other. Both mod- 
els are however much more efficient than the Smagorinsky model. In [14], we present a 
more complete study and show that same conclusions are obtained for the study of the 
evolution of the velocity relative error. For a smaller resolution 128X128, the anisotropic 
relaxation equations turn out to give better results than the hyperviscous parameterization. 

The previous computation is a typical merging of a small number of vortices. The 
initial configuration of these vortices leads to a rapid stabilization of a single structure. In 
order to study the effect of the parameterizations for a more complex dynamical evolution, 
we will consider an initial vorticity field composed by 50 vorticity patches, with random 
initial positions. We make the same study as in the previous case, comparing the efficiency 
of the different parameterizations in a low resolution computation (256X256), using as a 
reference a larger resolution computation (1024X1024). As in the previous case, the Navier- 
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Stokes or Smagorinsky parameterizations lead to a much less efficient computation than 
the hyperviscous or anisotropic relaxation equations ones. For instance, Fig. 7, for time 
T = 45, clearly shows that the Navier-Stokes equation do not even reproduce the qualitative 
properties of the vorticity fields for large time. This is due to the spurious effect of the 
viscous dissipation. The Smagorinsky model also leads to a too important diffusion. (As in 
previous computations, we alway use the smaller viscosity or parameter I, compatible with 
the stability of the computation. Numerical parameters are reproduced on table 4.) 

Fig. 8 shows the evolution of the velocity relative error for three low resolution (256X256) 
computations : hyperviscous (HV.Euler50.256), anisotropic relaxation (RELANI.Euler50.256) 
and anisotropic relaxation with local energy conservation (RELANILOC.Euler50.256). We 
use as a reference an hyperviscous computation with resolution 1024X1024. For smaller 
times, the three computations give equivalent results : a quasi-exponential evolution of the 
error, probably associated to the intrinsic instability of the initial condition. However, after 
a stage of several vortex merging, the anisotropic relaxation equations give much better 
results than the hyperviscous ones. This is illustrated by Fig.s 8 and 7. This last figure 
show that the position and the structures of the two final vortices is better reproduced. We 
propose the following interpretation : whereas the vortex motion is very sensitive to the 
initial condition, the result of the vortex merging depends less on the detail initial condi- 
tion. It rather depends on the dynamical invariants (the global vorticity, the energy). The 
relaxation equation which explicitly conserve the energy, and which describes the vorticity 
mixing in accordance with the statistical properties of the system lead to better results than 
the hyperviscous or the other parameterization studied here. We thus observe a statistical 
saturation of the error during this stage. 

These results show that the Navier-Stokes equation is very clearly the poorer model of 
two dimensional decaying turbulence. To compute the error, we have used, as a reference, 
an hyperviscous type computation. In order for the results to be non ambiguous, we have 
shown in [14] that all the conclusions we give here are also valid when a comparison is made 
with Navier-Stokes high resolution computations (HV.Euler50.1024) (see also figure 7). 

We conclude that in any situations, hyperviscous parameterizations and anisotropic re- 
laxation equations lead to much better results than Navier-Stokes equation or the Smagorin- 
sky model. The difference between hyperviscous and anisotropic relaxation equations com- 
putations is less important. However for small resolution (HV.Euler4.128) or for very com- 
plex flow evolution, the anisotropic relaxation equations give better results than the hyper- 
viscous parameterization. 

Our validation by comparison with higher resolution computation, does not allow us 
to compare different models for very large times. However, we think that the relaxation 
equation will then give much better results than the hyperviscous ones, because they ex- 
actly conserve the energy, because they treat vorticity mixing as far as possible according 
to the physical processes responsible for this mixing, and because they converge towards 
equilibrium states of the Euler equation. On the contrary hyperviscous model are not based 
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on any physical ground. 

We conclude by stating that similar results have been obtained in the context of the 
Quasi-Geostrophic dynamics (not reported, see [14]). Due to the slowing down of the dy- 
namical mixing of the potential vorticity for small Rossby deformation radius R (the inter- 
action is screened at a typical scale R), the relative efficiency of the anisotropic relaxation 
equations is even more evident in such a case. 

6 Discussion 

For geophysical flows, the scale at which the microscopic diffusion is relevant is so 
small that any numerical model of oceans or atmosphere implicitly assumes a small scale 
parameterization. In such physical situations, lots of phenomenon occur at the unresolved 
scales : convection, interaction with boundary, three dimensional turbulence at smaller 
scales, leading to molecular diffusion, etc ... For this reason, the fundamental equation to 
start with to describe such flows is not even clear [1]. It is however generally recognized 
that the dynamics of geophysical flows is two-dimensional like (inverse energy cascade) and 
that the dynamics of the potential vorticity plays a major role. In order to parameterize 
the potential vorticity mixing, we have considered the simplest model having the main 
properties of large scale geophysical flows : the barotropic Quasi-geostrophic model. As the 
equations are formally equivalent, we have led our study mainly in the 2D incompressible 
Euler-equation framework. All the results obtained are valid for both of these models, and 
may be straightforwardly generalized to the multi-layered Quasi-Geostrophic model. The 
problem we have addressed is to describe the evolution of the coarse-grained vorticity, the 
fundamental equation being the Euler (or QG) equations. In this study, we have considered 
a regime of freely decaying turbulence, with vortex merging and filamentation. As we mainly 
study this transient behavior, this study also apply for the two-dimensional Navier-Stokes 
equations (or viscous QG) with a molecular viscosity acting on a scale much smaller than 
the resolved one. 

In section 2, we have first defined the coarse-graining. We have then shown that the 
parameterization problem is equivalent to the determination of a turbulent current. We 
have obtained at leading order a deterministic approximation of the turbulent current : the 
action of the strain tensor on the gradient of the vorticity[16, 15, 14]. We have shown using 
numerically obtained vorticity fields that this approximation is a good one in the regime 
of vortex merging and filamentation. The leading contribution comes from the part of the 
turbulent current which is absent in the case of an ensemble average (like in usual Reynolds 
tensor). This shows that the turbulent current is dominated by systematic correlation be- 
tween resolved and unresolved scale. This essentially deterministic process is qualitatively 
associated to structures (filaments) passing from resolved to unresolved scales). This deter- 
ministic approximation exactly conserves the energy. 

This deterministic approximation acts on the vorticity field as a positive viscosity in 
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the stretching direction of the strain tensor, and as a negative viscosity in the expansion 
direction. Due to this last contribution, a numerical model using directly this approxima- 
tion would lead to numerical instabilities. In order to obtain stable equations, respecting 
the actual diffusion in the contraction direction, and the energy conservation, we look for 
equations maximizing the mixing entropy for the vorticity, while imposing these dynamical 
constraints. We obtain anisotropic relaxations equations which have the further property 
to relax towards the statistical equilibrium of the Euler-equation. 

Some entropy based parameterization of the turbulent current have been proposed pre- 
viously in the case of topography dominated flows, by Holloway (see [22] for a review or 
[23] for some application in the context of oceanographic flows). The maximum entropy 
production principle has been proposed in [12], but an isotropic turbulent current was then 
imposed. The maximum entropy production principle has no theoretical or physical ground. 
The anisotropic relaxation equations do not for instance describe correctly what happens in 
the strain tensor stretching direction. We use here this principle as a trick to obtain stable 
equations, with given constraints, and respecting the statistical tendency of the vorticity 
mixing. 

In order to show the interest of the anisotropic relaxation equations, we have compared 
its efficiency with the Navier-Stokes equations, or other usual parameterizations such as the 
Smagorinsky model or the hyperviscous model. We have compared the results with higher 
resolution Navier-Stokes or hyperviscous computations. As they overestimate strongly the 
diffusion, the Navier-Stokes or the Smagorinsky model give clearly the worst results. The 
anisotropic relaxation equations give better results than the hyperviscous model, for com- 
plex flow evolution, for smaller resolutions, or for the Quasi-Geostrophic model. 

The results we have obtained clearly show that the anisotropic relaxation equations 
allows an increase of the predictability with respect to other parameterizations. We note 
however that we have suppose a complete knowledge of the initial condition. The real prob- 
lem of predictability for geophysical flows is much more complex, as lots of error sources 
have to be taken into account. We have studied only the effect of the parameterization of 
the potential vorticity mixing. 

Let us conclude on a discussion on the limitations of this work and on its possible 
interest in a wider context. We have First we have studied the regime of vortex merging 
and strong filamentation. However, for inviscid two-dimensional freely decaying turbulence, 
on the latter stage of the evolution of the flow, the coarse-grained field is close to a stationary 
state. In this different regime, no more structures are transported from larger to smaller 
scales (scale separation). In such a case, the main contribution to the turbulent current 
probably comes from the stochastic part of the current. The effect of these fluctuations 
is a very slow modification of the structure of the quasi-stationary coarse-grained flow. 
An ensemble average is probably the best way to address such a relaxation. We note the 
application of a quasi-linear theory in [11] which may be of interest for this regime. 

In the situation we have studied, because the processes are very rapid, the effect of 
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a very small viscosity (acting on scales much smaller than the resolution one) is actually 
negligible. However, in the stage of a quasi-stationary state relaxation, the effect of a very 
small viscosity would probably be very important, as it should determine the distribution 
of the unresolved scale vorticity. A statistical mechanics based approach is proposed in 
[26]. There a crude closure, using the simplest energy-enstrophy statistical theory for flow 
with topography, is proposed. In [25] the effect of a small forcing and dissipation is also 
considered. 

Another improvement of our work may be done by the description, not only of the av- 
eraged but also of the higher vorticity moments. In [24] an equation hierarchy is derived, to 
describe these moment evolution. This hierarchy is then closed using a maximum entropy 
production principle. Such ideas have been applied in the context of a barotropic ocean 
model [25]. 

The current computational power allows to make very high resolution simulations of 
the two-dimensional Navier Stokes equations. In such a case, the precise simulation of rapid 
physical phenomena does not necessitate any parameterization. However, for more complex 
physical situations, for very long time simulations or for more complex geophysical flow 
models, the necessity to use good small scales parameterizations is evident. Let us give 
some examples. In the modeling of Jupiter's troposphere, we have shown that the very 
small value of the Rossby deformation radius renders the potential vorticity dynamics very 
slow [27]. The very long time then needed, for the flow to organize, does not allow to obtain 
the characteristic very strong jets typical of Jupiter's vortices, when using Navier-Stokes 
equations. On the contrary, using relaxation equations permits to model these features 
very precisely, even using quite low resolution models [27]. For ocean modeling, let us 
discuss the example of the Zapiola anticyclone, in the Argentina basin, in South Atlantic. 
In the late 90's, Topeix-Poseidon data has clearly revealed this very strong anticyclonic 
structure, whereas up to date South Atlantic modeling where not able to reproduce it, or 
were reproducing a structure much weaker than the observed ones [28]. In [28], for instance, 
the anticyclone modeling is convincing, using the up to date model resolution. Moreover, 
for global ocean models, or/and when much longer computational times are needed, for 
instance for climate modeling, the correct description of such structures would necessitate 
unavailable computational power. Moreover, as shown by Fig. 7 of this study, for very 
long time scales, the use of rough parameterization alters the fluctuations of vorticity, or 
potential vorticity. 

The modeling of geophysical flows actually require the best possible parameterization. 
The analysis we have done could be straightforwardly generalized to the case of multi- 
layered Quasi-geostrophic models. For z-coordinates models, or isopicnal models in use in 
geophysical flow modeling, the generalization of the ideas developed in this article should be 
considered ; the problem being much more simple in models clearly identifying the potential 
vorticity as a key dynamical variable. We conclude by saying that the modeling of the poten- 
tial vorticity mixing in geophysical fluid dynamics is only one part of the parameterization 
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problem, but probably one of the more essential one. 
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Tab. 1 - Comparison of the contribution of the divergence of the components of the tur- 
bulent current, evaluated with several norms ||/|| 2 = ^/l/ \ D\ J D dr f 2 , IHI^ = Sup D {|/|}, 
H/l^ = 1/ \D\ J D dr \ f\. The divergence of the turbulent current V.J. = V.J, + V.J S 
is dominated by the contribution of its deterministic part V.J^. The stochastic part 
V.J S = V. (<^u) + V. ^liuj is dominated by V. (^G) . This last term would be the only 
contribution in the case of an ensemble average. 
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Tab. 2 - Norm of the difference between the divergence of the deterministic part of the 
turbulent current V.J, and its first and second order approximations. The leading order 
approximation of V.J, is correct up to an error of order 10%, for the three norms considered. 
This error is of the same order as the one corresponding to the stochastic part V. J s 
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FlG. 1 - Upper picture : a typical vorticity field uj for freely decaying turbulence in the 
regime of filamentation and vortex coalescence (resolution 1024X1024). Lower picture : 
the corresponding averaged vorticity field u, for a Gaussian homogenization at the scale 
I = 27r/128 ~ 0.049. Fig. 2 shows the corresponding turbulent current. 
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T = 14.00 




Fig. 2 - The divergence of the turbulent current V.J W , computed for a Gaussian homog- 
enization at a scale I = 27r/128 ~ 0.049, for the vorticity field shown on Fig. 1. The 
comparison with Fig. 1 illustrates that important values of the divergence of the turbu- 
lent current are associated with structures with typical scale I, especially filaments. Not all 
structures at scale I give a strong contribution, but only the ones correlated with a strong 
mean strain (see Fig. 3). 
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Fig. 3 - The mean strain a (positive eigenvalue of the symmetric part of the strain tensor) 
for the velocity field corresponding to the homogenized vorticity uJ (Fig. 1). We note the 
correlation between areas of strong strain, areas with strong vorticity gradients (Fig. 1) and 
areas where the divergence of the turbulent current is high (Fig. 2). 



Calcul 


Resolution 


V ou I 


dt 


Re 


DNS.Euler4.1024 


1024x1024 


1.57 10" 5 


6.14 10" 4 


400 000 


REL.Euler4.256 


256x256 


6.28 10" 4 


2.45 10" 3 


10 000 


RELVAR.Euler4.256 


256x256 


1.26 


2.45 10" 3 




RELANI.Euler4.256 


256x256 


1.33 10~ 2 


2.45 10" 3 




RELANILOC.Euler4.256 


256x256 


1.33 10" 2 


2.45 10" 3 




DNS.Euler4.256 


256x256 


6.28 10" 4 


2.45 10" 3 


10 000 


SMA.Euler4.256 


256x256 


1.33 10" 2 


2.45 10" 3 




HV.Euler4.256 


256x256 


1.92 10" 8 


2.45 10" 3 




RELANI.Euler4.128 


128x128 


2.45 10" 2 


4.91 10" 3 




HV.Euler4.128 


128x128 


1.53 10 -7 


4.91 10" 3 





Tab. 3 - Numerical parameters for the computations of the coalescence of four vortices 
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Fig. 4 - Vorticity evolution for the coalescence of four vortices. This computation uses a 
Navier-Stokes model, with resolution 1024x1024 (DNS.Euler4.1024). Time goes from T = 
to T = 99. We consider this computation as representative of the inertial limit (energy loss 
of 0.7% between T = and T = 100). 
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Fig. 5 - Velocity relative error on the low resolution computation, for the coalescence of 
four vortices, for parameterizations of small scales with respectively : a constant diffusivity 
(Navier-Stokes) (DNS.Euler4.256), isotropic relaxation equations with constant diffusion 
(REL.Euler4.256), anisotropic relaxation equations (RELANI.Euler4.256) and anisotropic 
relaxation equations with local energy conservation (RELANILOC.Euler4.256). 



Calcul 


Resolution 


V ou I 


dt 


Re 


DNS.Euler50.1024 


1024x1024 


3.14 10" 5 


6.14 10" 4 


200 000 


DNS.Euler50.1024 


1024x1024 


5.99 10" 11 


6.14 10" 4 




RELANI.Euler50.256 


256x256 


1.33 10" 2 


2.45 IQ-' 6 




RELANILOC.Euler50.256 


256x256 


1.33 10" 2 


2.45 10" a 




HV.Euler50.256 


256x256 


1.92 10" s 


2.45 10" 3 





Tab. 4 - Numerical parameters for the various parameterization of the coalescence of 50 
vorticity patches. 
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Erreur relative sur la vorticit 



36e-3 
32e-3 - 
28e-3 
24e-3 
20e-3 - 
16e-3 
12e-3 

8e-3 -| 

4e-3 




Er_w 



RELAN1LOC.EULER4.256 

SMA.EULER4.256 

HV.EULER4.256 




Fig. 6 - Same as Fig. 5, but for the Smagorinsky model (SMA.Euler4.256), the 
anisotropic relaxation equation (RELANI.Euler4.256) and an hyperviscous parameteriza- 
tion (HV.Euler.4.256) ; for computations with resolution 256x256. 
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Fig. 7 - Comparison of vorticity field, obtained at time T = 45.01, from the same initial 
condition, for two different high resolution computations : hyperviscous (HV.Euler50.1024) 
and Navier Stokes (DNS. Euler50. 1024) ; and four low resolution computations with differ- 
ent small scale turbulence parameterization : hyperviscosity (HV.Euler50.256), anisotropic 
relaxation equations (RELANI.Euler50.256), anisotropic relaxation equations with local 
enrgy conservation (RELANILOC.Euler50.2sB) and Navier Stokes (DNS.Euler50.256). 




Fig. 8 - Velocity relative error versus time, for three low resolution computations (256x256), 
with different small scale turbulence parameterization : anisotropic relaxation equations 
(RELANI.Euler50.256), hyperviscous model (HV.Euler50.256) and anisotropic relaxation 
equations with local energy conservation (RELANILOC.Euler50.256). The initial condition 
is made of 50 vorticity patches with random positions. The reference computation is a 
1024X1024 resolution hyperviscous one (HV.EULER50.1024) 
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